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Air shower simulations are essential for interpreting data from cosmic ray experiments. At highest energies 
though, a microscopic treatment of a whole shower is not possible any more, since it would require a huge amount 
of CPU-time. We review hybrid approaches of air shower simulation which try to overcome this problem without 
giving rise to artificial fluctuations as generated by the thinning algorithm. 



1. Introduction 

Ultra-high energy cosmic rays (UHECR) are 
currently of great interest. One still ignores the 
nature of these particles, acceleration mechanisms 
and possible sources. Direct measurements at 
these energies are impossible due to a very low 
flux. Upon entering the atmosphere, cosmic rays 
induce extensive air showers (EAS), cascades of 
particles produced by collisions of the primary 
and subsequent secondaries with air nuclei. Ex- 
periments measure these showers and try to de- 
duce properties of the primary particle from prop- 
erties of the shower. One can therefore consider 
the atmosphere as a huge detector in which a pri- 
mary cosmic ray is absorbed and observed. 

There are basically two types of experimental 
set-ups. Ground arrays (i.e. AGASA, Auger) 
measure the density of charged particles on the 
ground. The lateral distribution function (LDF) 
is studied to deduce the energy. One can estimate 
the arrival direction from the different trigger- 
times of the detectors as the shower front passes 
through the array. The other kind are optical 
detectors. They collect the light emitted by fluo- 
rescence of nitrogen molecules, which are excited 
by the flux of charged particles going through 
the atmosphere. They measure therefore directly 
the longitudinal profile of a shower. Both types 
of experiments have advantages and drawbacks. 
Ground arrays have a high duty cycle and are not 
very sensitive to details of the atmosphere. Flu- 
orescence experiments depend somewhat less on 



models, but are influenced by atmospheric fluc- 
tuations and have a duty cycle of typically 10%. 

For both set-ups, reconstruction of primary 
properties depends on how good one understands 
the interactions in the atmosphere. Air shower 
simulations are therefore crucial for cosmic ray 
physics as they are needed for interpretation of 
the data. A major difficulty is that no experi- 
mental data from accelerators is available at these 
energies. LHC is still three orders of magni- 
tudes (lab. system) lower than the currently high- 
est energies measured. Another problem is that 
air shower development is dominated by forward 
scattering, since most energy is carried by the 
leading particles, but accelerators measure mostly 
at mid-rapidity. Also, the targets are light nuclei, 
commonly less well studied. Therefore, physics 
has to be extrapolated into unknown regions (en- 
ergy and phase-space) which makes interpretation 
of data less precise. On the other hand, one can 
argue that cosmic rays provide a unique oppor- 
tunity to study physics at high energies that may 
never be reached by accelerator experiments. 

Air showers develop rapidly in the atmosphere. 
A rule of thumb is that the number of charged 
particles at the maximum of the longitudinal pro- 
file is about 60% of the energy measured in GeV; 
a 10 11 GeV shower has about 60 billion parti- 
cles. Hence, microscopic simulations seem quite 
impossible. E.g., a 10 10 GeV shower would take 
about a year to compute on current CPUs. The 
thinning algorithm pQ reduces the CPU time but 
introduces artificial fluctuations. Hybrid Simula- 
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tions 2 try to solve this problem by following 
the high energy part of an air shower in detail, 
and using efficient approximations below a given 
threshold. 

2. The thinning problem in microscopic 
treatment of air shower simulations 

Current air shower simulations use the thin- 
ning algorithm in order to limit the computa- 
tion time. The original idea is to follow only one 
particle of all secondaries below a given thresh- 
old E t h = fthEo, chosen with the probability 
Pi = Ei/Etot- Neglected particles are compen- 
sated by attributing a higher weight Wi = 1/pi 
to the chosen particle. In principle, the probabil- 
ity of some particle below the threshold to be fol- 
lowed can be an arbitrary function, but the weight 
attributed has to be 1/pi- Statistical thinning 
means pi — Ei /E to t (as above, but the number of 
followed particles is not fixed to one) . Further re- 
finements are possible by imposing a weight limit 
on the algorithm, or not to thin beyond a given 
distance from the shower axis. The disadvantage 
of thinning is that it introduces artificial fluctua- 
tions, and one has to be careful to control these. 

The hybrid approach tries to overcome these 
problems by computing the sub-showers of parti- 
cles efficiently instead of discarding these. 

3. Hybrid approaches 

As mentioned, hybrid methods replace sub- 
threshold particles with the sub-showers induced 
by them. We are going to discuss two different 
approaches, where these sub-showers are obtained 
from a shower library or the solution of cascade 
equations. They differ also in the fact, that the 
former replaces a fluctuating shower whereas the 
latter solves for mean sub-showers. 

3.1. Bartol approach to hybrid simulations 

An obvious choice is to implement a data base 
of pre-simulated showers. During the simulation, 
all sub-showers below a given energy threshold are 
replaced with a sample from this data base. The 
Bartol [Jj approach consists of a shower library 
of pion-initiated showers for different injection 
depths, energies and zenith angles. Nucleons are 
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Figure 1. Principle of the hybrid approach: com- 
pute sub-showers with efficient algorithms. 



followed by Monte-Carlo method, whereas kaons 
are treated to behave as pions with respect to in- 
teraction. The resulting longitudinal profiles are 
then fitted to a Gaisser-Hillas function, and the 
parameters are recorded. In addition, the num- 
ber of muons at ground level for different energy 
thresholds is stored. The library itself is built 
with a boot-strap method: Starting from low en- 
ergy, the pre-simulated showers are used for the 
computation of higher energies. During the sim- 
ulation of a given shower, each particle falling be- 
low the threshold E t h = O.OIE'o is replaced with 
a sub-shower sampled from the data base. Note 
that these sub-showers are not average profiles 
but account fully for fluctuations even below the 
hybrid threshold. 
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3.2. Hybrid approach with cascade equa- 
tions 

Cascade equations are one-dimensional trans- 
port equations which are solved numerically. 
The solution is a mean shower profile without 
any fluctuations. Given a number of particles 
h n (E, X)dE, of type n, at given altitude X and 
energy E, the probability for interaction and de- 
cay is 

with A being the mean free path of the particle as 
a function of the energy, p being the density of the 
air, and d n = m n /(cT n ) being the decay constant 
(E/d n is the decay length in the lab. system). 

Accounting for particles produced at higher 
energies gives rise to the following system of 
hadronic cascade equations 5 : 
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The functions W mn (E' ', E) are the energy-spectra 
j£ of secondary particles of type n in a colli- 
sion of hadron m with air-molcculcs. D rnn (E' ', E) 
are the corresponding decay-functions. Equation 
(|2j is a typical transport equation with a source 
term. The first term with the minus-sign accounts 
for particles disappearing by collisions or decays, 
whereas the source term accounts for production 
of secondary particles by collisions or decays of 
particles at higher energies. 

The initial condition for the primary cosmic ray 
is given by 



h n {E, X — X. m ) — 6 nm S(E — E m ). 
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When using the hybrid approach, the initial con- 
dition is a superposition of all particles E < E t h 
produced above the threshold. 

3.3. Low energy source functions 

At low energies, the three-dimensional spread 
of particles becomes important. For the computa- 
tion of the longitudinal profile only, one can apply 



corrections to the one-dimensional cascade equa- 
tions in order to account for neglecting the lateral 
expansion. This is done in the CONEX model @]. 

For the computation of lateral distribution 
functions one can switch back to the Monte-Carlo 
scheme. Particles are generated from the so- 
called source function. 
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Source functions are used in the SENECA model 
|B] and by Dedenko et al. |3j . Particles are gener- 
ated according to Q and placed along the shower 
axis. 

The source function contains a certain energy 
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which should be used to produce secondary parti- 
cles. By choosing a fraction / of this energy, one 
determines the weight of the secondaries. Ideally 
this should be equal to one, but this would still 
take too much of CPU power. However, once de- 
termined, the weight stays constant even in the 
subsequent tracking with MC method. This is an 
advantage over the thinning method where the 
weight of particles is more difficult to control. 

The choice of the energy threshold E max , where 
to switch back from CE to MC is crucial. A lower 
value is desirable for computation speed, whereas 
a higher value might be needed to achieve the 
desired precision for the lateral distribution func- 
tion. The best value to choose depends therefore 
also on the observable to be computed. For ex- 
ample, for longitudinal profiles the lateral spread 
of muons is not important, and one can choose a 
lower transition energy. 

3.4. Electromagnetic showers 

The electromagnetic part of the shower can 
be treated in a similar way. In SENECA, pre- 
simulated sub-showers for a 2.5 g/cm 2 thick layer 
of air are stored in a table. The energies are dis- 
cretized in ten bins per decade E t = 10 1 / 10 . The 
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Figure 2. A single mean CE shower compared 
with 1000 averaged MC-showers (10~ 5 thinning). 



table is then a matrix V™ n , where the indices i,j 
represent primary and secondary energies, m, n 
stand for the particle types, photons and elec- 
trons/positrons. If we have g"(X) particles of 
type n and energy Ei, the corresponding spec- 
trum at X + AX is 

gUX + AX) = Y,V™97(X) • 

Dedenko et al., Lagutin et al. [5], as well as 
CONEX choose to implement cascade equations in 
a similar way to (J2J) for all electromagnetic in- 
teractions, i.e. pair production, bremsstrahlung, 
etc. At low energies, CONEX implements a higher 
effective path length, in order to correct for the 
neglected transverse dimensions. In SENECA one 
can apply a table of pre-simulated sub-showers. 
This way, the longitudinal profile can be calcu- 
lated in a quick way, using the slow MC method 
only for initial fluctuations. 

4. Results and Comparisons 

4.1. Internal checks 

Precision is crucial for the hybrid approach. 
Any lack of precision will appear in the solution 
systematically. Therefore it is necessary to con- 
firm whether CE and MC give the same results. 
Since the cascade equations give as result a mean 



Figure 3. Lateral distributions for elec- 
trons/positrons, muons and photons. The cut-off 
energies are 50 MeV for muons and 1 MeV for the 
EM-particles. 



shower, we compare the longitudinal profile of one 
CE shower to the average of 1000 MC generated 
showers. Here we do not use MC for the initial 
part, i.e. / = 1. The agreement of the two curves 
is shown in Fig. |21 

Since the CE approach is ideal to compute an 
average shower, one might ask how good the max- 
imum of a mean shower describes the mean of 
fluctuating showers, as measured by experiments. 
Assuming the shower profiles are described by 
functions fi(x), the mean shower is 

1 N 

»=i 

The maximum is defined by f'(x) — so 
f'(X max ) = and //(X maX)i ) = 0. Expanding 
each profile i around its maximum yields 

fi\*E) — /i(-^max,i) 

. fj (Xmaxj) , , 2 

+ ... (6) 

Differentiating to find the maximum gives 

= f'(X max ) = — y^//(A max ) 
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Figure 4. Energy spectra for both simulation 
methods. 



jy ^ ] fi (^max.ijt^max X maK i) ) (7) 

which leads to 

X max = (X maXt i) for /"(X maXj i) = c . (8) 

So, if the curvature at X max i is approximately 
a constant, the maximum of a mean shower cor- 
responds to the mean of maxima. Simulations 
show that this relationship JHJ) is valid within a 
few g/cm 2 . Hence, for a quick estimate one can 
use a mean shower computed with the CE ap- 
proach, instead of doing hundreds of fluctuating 
showers. 

The same does not hold for the shower size, 
S'max- Evaluating formula © at X maK yields 

+ ^(^L,,)-^,,) 2 ) (9) 

S max of a mean shower is smaller than the mean 
of S^nax of many showers, since c is negative. Be- 
cause of energy conservation, we expect the mean 
shower to be somewhat wider than on average 
(larger A in the Gaisser-Hillas formula). In that 
sense, a mean shower is not a typical shower, but 
it can be used to get a good estimate of the mean 

Xma,x- 
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Figure 5. The spread of X max as a function of the 
hybrid threshold /. 



4.2. Lateral Distribution Functions 

A check of lateral distribution functions (LDF) 
is shown in Fig. The densities for elec- 
trons/positrons, muons, and photons agree with 
the MC simulated spectra of a vertical 10 19 eV 
proton induced shower. In this example we com- 
pare two fluctuating showers as opposed to com- 
paring the average behavior. The high energy 
part of the shower is calculated in the same way 
in both the MC and the CE simulation. Tech- 
nically this is achieved by implementing a high 
energy particle stack and using the same seed for 
the pseudo-random number generator. The high 
energy stack takes care that the sequence of ran- 
dom numbers is exactly the same until the first 
particle below the threshold appears in the calcu- 
lation. 

An important parameter for the LDFs is the 
threshold in energy for switching back to the MC 
method. We choose 10 GeV for electromagnetic 
cascades and 10 TeV for the hadronic part. This 
relatively high value is necessary to reproduce 
correctly the lateral spread of the muons. 

4.3. Energy Spectra 

Energy spectra of all particles are directly cal- 
culated in the CE. An example for energy spectra 
at 600 g/cm 2 is shown in Fig. 0] Note the dif- 
ferences in K-short and K-long below 10 5 GeV; 
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Figure 6. Mean A max of conex and 
SENECA compared with CORSIKA results. The 
Seneca proton result for NeXus 3.97 has been 
computed with CE only ( A max of a mean shower) . 



this is the energy region where decays of K-shorts 
start dominate over interaction. 

4.4. Fluctuations 

An important check for the hybrid approach is 
whether it can reproduce natural fluctuations in 
the A max distribution. Fig. [S] shows the spread 
{ ) 2 as a function of the 



threshold / (below E = JEq cascade equations 
are used). As of / = 0.01 the spread seems to con- 
verge. Interestingly, already / = 0.99 reproduces 
90% of the spread, i.e. the very first interaction 
is responsible for most of the fluctuations. 

4.5. Shower maximum 

A comparison with the well tested CORSIKA [5] 
model is given in Fig. Since the physics 
content in terms of external models is the 
same, the hybrid approaches should give simi- 
lar results. QGSJET01 [HI], SIBYLL 2.l[TI] and 
neXus 3.97^21 are use d as high energy hadronic 
models in the framework of conex and SENECA. 

When computing an iron induced shower with 
CE, one does not need any additional tables for 
nucleus primaries, if the energy per nucleon is 
above the initial fluctuation threshold. For di- 
rect computation of a mean shower, one would 
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Figure 7. Arrival time distributions of photons 
and the corresponding maximum weight in a 
given bin. The peak in the MC figure is due to a 
single particle with a huge weight. 



have in principle to do tables for all possible nu- 
cleus primaries up to A = 56. In that case, it 
is more reasonable to average over many iron-air 
collisions (such that no nuclei are left) before solv- 
ing the transport equations. 

4.6. Arrival time 

Arrival times are interesting, since modern ex- 
periments try to extract information about the 
primary via e.g. the rise time. Here, the thin- 
ning procedure has a great disadvantage, since 
the weight of a given particle can be dominant for 
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Figure 8. A 80° inclined 10 19 eV proton induced 
shower for pure MC and the hybrid (CE) method. 



the signal, as shown in Fig. [7] The largest peak is 
due to a single particle with a huge weight. Par- 
ticles generated from the source function have a 
constant weight, which is adjustable. However, 
refined thinning algorithms have a weight limit, 
with which the maximum weight of a particle can 
be defined. 

This is one example where the technique of CE 
has a useful side-effect, other than the main ad- 
vantage of reduced CPU-time. 

4.7. Very inclined showers 

Inclined showers are of great interest for wa- 
ter tank detectors as in the Auger observatory. 
These are very different from near-vertical show- 
ers, due to a huge path-length. For 80° inclined 
showers, we have 6 times the thickness of the at- 
mosphere, and at 90° 36 times. The electromag- 
netic part is almost absorbed in the atmosphere 
after some 2000 g/cm 2 ; only muons continue to 
propagate with few interactions due to energy- 
loss, decays, bremsstrahlung, etc. The result is a 
relatively flat profile, see Fig. |SJ The accompa- 
nying electrons/positrons and photons come from 
interactions and decays of muons. When calcu- 
lating such a profile with CE, we solve only up to 
2000 g/cm 2 slant depth, since beyond that value 
muons dominate the shower. The profiles from 
the MC and CE approaches agree nicely. But 



since the distance is large, a small error in the cre- 
ation of hadrons from the source function could 
result in a wrong lateral shape. This is shown in 
Fig. [5J where we plot the muon density at ground 
level in the shower plane. The deflection of posi- 
tive and negative muons due to the geomagnetic 
field is clearly visible. This is important to take 
into account when analyzing experimental data 
from ground arrays. Fig. lUtop) shows the MC 
event and Fig. [5{middle) a corresponding result 
with CE. Since the particles are placed with zero 
angle along the shower axis, the pattern looks dis- 
tinctly different. In Fig. ^bottom) a hybrid cal- 
culation is shown, where a mean transverse mo- 
mentum is applied to all particles generated by 
the source function. A value of pt ~ 0.3 GeV 
seems to be sufficient to reproduce the right pat- 
tern. 

5. CPU times 

In ref. [H] it was shown that the CE approach is 
a factor 20-40 faster than the MC approach, when 
asking for LDFs of the same statistical quality. A 
similar enhancement was found in ref. 7; by using 
the shower library approach. 

Of course, it depends very much on the observ- 
able to be computed when comparing CPU times. 
But in general, a typical longitudinal profile with- 
out any lateral spread can be computed in about 
a minute on a lGhz CPU. A useful lateral distri- 
bution function takes about 10 minutes. 

6. Conclusions 

Hybrid calculations allow one to reduce con- 
siderably the computation time for air shower 
simulations. The natural fluctuations arise from 
the first few interactions in the atmosphere and 
are computed in traditional Monte-Carlo method. 
Below a given threshold, particles are then re- 
placed with sub-showers taken from a shower li- 
brary or by the solution of cascade equations. 
The CE can be solved down to lowest energies, 
if corrections are applied to account for the ne- 
glected lateral expansion. The precise lateral 
spread can be computed again by MC method, 
with low energy particles created from the source 
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function. 

The hybrid approach, i.e. the combination of 
traditional Monte-Carlo with efficient numerical 
methods provides a powerful tool for studying 
ultra-high energy cosmic rays. 
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Figure 9. The density of muons in the shower 
plane. Top panel: MC result, middle panel: CE 
result, lower panel: CE result with p t -kick in 



